home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / gnu / emacs.lha / emacs-19.16 / lisp / float.el < prev    next >
Lisp/Scheme  |  1993-06-09  |  15KB  |  458 lines
  1. ;;; float.el --- floating point arithmetic package.
  2.  
  3. ;; Copyright (C) 1986 Free Software Foundation, Inc.
  4.  
  5. ;; Author: Bill Rosenblatt
  6. ;; Maintainer: FSF
  7. ;; Keywords: extensions
  8.  
  9. ;; This file is part of GNU Emacs.
  10.  
  11. ;; GNU Emacs is free software; you can redistribute it and/or modify
  12. ;; it under the terms of the GNU General Public License as published by
  13. ;; the Free Software Foundation; either version 2, or (at your option)
  14. ;; any later version.
  15.  
  16. ;; GNU Emacs is distributed in the hope that it will be useful,
  17. ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
  18. ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19. ;; GNU General Public License for more details.
  20.  
  21. ;; You should have received a copy of the GNU General Public License
  22. ;; along with GNU Emacs; see the file COPYING.  If not, write to
  23. ;; the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  24.  
  25. ;;; Commentary:
  26.  
  27. ;; Floating point numbers are represented by dot-pairs (mant . exp)
  28. ;; where mant is the 24-bit signed integral mantissa and exp is the
  29. ;; base 2 exponent.
  30. ;;
  31. ;; Emacs LISP supports a 24-bit signed integer data type, which has a
  32. ;; range of -(2**23) to +(2**23)-1, or -8388608 to 8388607 decimal.
  33. ;; This gives six significant decimal digit accuracy.  Exponents can
  34. ;; be anything in the range -(2**23) to +(2**23)-1.
  35. ;;
  36. ;; User interface:
  37. ;; function f converts from integer to floating point
  38. ;; function string-to-float converts from string to floating point
  39. ;; function fint converts a floating point to integer (with truncation)
  40. ;; function float-to-string converts from floating point to string
  41. ;;                   
  42. ;; Caveats:
  43. ;; -  Exponents outside of the range of +/-100 or so will cause certain 
  44. ;;    functions (especially conversion routines) to take forever.
  45. ;; -  Very little checking is done for fixed point overflow/underflow.
  46. ;; -  No checking is done for over/underflow of the exponent
  47. ;;    (hardly necessary when exponent can be 2**23).
  48. ;; 
  49. ;;
  50. ;; Bill Rosenblatt
  51. ;; June 20, 1986
  52. ;;
  53.  
  54. ;;; Code:
  55.  
  56. ;; fundamental implementation constants
  57. (defconst exp-base 2
  58.   "Base of exponent in this floating point representation.")
  59.  
  60. (defconst mantissa-bits 24
  61.   "Number of significant bits in this floating point representation.")
  62.  
  63. (defconst decimal-digits 6
  64.   "Number of decimal digits expected to be accurate.")
  65.  
  66. (defconst expt-digits 2
  67.   "Maximum permitted digits in a scientific notation exponent.")
  68.  
  69. ;; other constants
  70. (defconst maxbit (1- mantissa-bits)
  71.   "Number of highest bit")
  72.  
  73. (defconst mantissa-maxval (1- (ash 1 maxbit))
  74.   "Maximum permissible value of mantissa")
  75.  
  76. (defconst mantissa-minval (ash 1 maxbit)
  77.   "Minimum permissible value of mantissa")
  78.  
  79. (defconst floating-point-regexp
  80.   "^[ \t]*\\(-?\\)\\([0-9]*\\)\
  81. \\(\\.\\([0-9]*\\)\\|\\)\
  82. \\(\\(\\([Ee]\\)\\(-?\\)\\([0-9][0-9]*\\)\\)\\|\\)[ \t]*$"
  83.   "Regular expression to match floating point numbers.  Extract matches:
  84. 1 - minus sign
  85. 2 - integer part
  86. 4 - fractional part
  87. 8 - minus sign for power of ten
  88. 9 - power of ten
  89. ")
  90.  
  91. (defconst high-bit-mask (ash 1 maxbit)
  92.   "Masks all bits except the high-order (sign) bit.")
  93.  
  94. (defconst second-bit-mask (ash 1 (1- maxbit))
  95.   "Masks all bits except the highest-order magnitude bit")
  96.  
  97. ;; various useful floating point constants
  98. (setq _f0 '(0 . 1))
  99.  
  100. (setq _f1/2 '(4194304 . -23))
  101.  
  102. (setq _f1 '(4194304 . -22))
  103.  
  104. (setq _f10 '(5242880 . -19))
  105.  
  106. ;; support for decimal conversion routines
  107. (setq powers-of-10 (make-vector (1+ decimal-digits) _f1))
  108. (aset powers-of-10 1 _f10)
  109. (aset powers-of-10 2 '(6553600 . -16))
  110. (aset powers-of-10 3 '(8192000 . -13))
  111. (aset powers-of-10 4 '(5120000 . -9))
  112. (aset powers-of-10 5 '(6400000 . -6))
  113. (aset powers-of-10 6 '(8000000 . -3))
  114.  
  115. (setq all-decimal-digs-minval (aref powers-of-10 (1- decimal-digits))
  116.       highest-power-of-10 (aref powers-of-10 decimal-digits))
  117.  
  118. (defun fashl (fnum)            ; floating-point arithmetic shift left
  119.   (cons (ash (car fnum) 1) (1- (cdr fnum))))
  120.  
  121. (defun fashr (fnum)            ; floating point arithmetic shift right
  122.   (cons (ash (car fnum) -1) (1+ (cdr fnum))))
  123.  
  124. (defun normalize (fnum)
  125.   (if (> (car fnum) 0)            ; make sure next-to-highest bit is set
  126.       (while (zerop (logand (car fnum) second-bit-mask))
  127.     (setq fnum (fashl fnum)))
  128.     (if (< (car fnum) 0)        ; make sure highest bit is set
  129.     (while (zerop (logand (car fnum) high-bit-mask))
  130.       (setq fnum (fashl fnum)))
  131.       (setq fnum _f0)))            ; "standard 0"
  132.   fnum)
  133.       
  134. (defun abs (n)                ; integer absolute value
  135.   (if (>= n 0) n (- n)))
  136.  
  137. (defun fabs (fnum)            ; re-normalize after taking abs value
  138.   (normalize (cons (abs (car fnum)) (cdr fnum))))
  139.  
  140. (defun xor (a b)            ; logical exclusive or
  141.   (and (or a b) (not (and a b))))
  142.  
  143. (defun same-sign (a b)            ; two f-p numbers have same sign?
  144.   (not (xor (natnump (car a)) (natnump (car b)))))
  145.  
  146. (defun extract-match (str i)        ; used after string-match
  147.   (condition-case ()
  148.       (substring str (match-beginning i) (match-end i))
  149.     (error "")))
  150.  
  151. ;; support for the multiplication function
  152. (setq halfword-bits (/ mantissa-bits 2)    ; bits in a halfword
  153.       masklo (1- (ash 1 halfword-bits)) ; isolate the lower halfword
  154.       maskhi (lognot masklo)        ; isolate the upper halfword
  155.       round-limit (ash 1 (/ halfword-bits 2)))
  156.  
  157. (defun hihalf (n)            ; return high halfword, shifted down
  158.   (ash (logand n maskhi) (- halfword-bits)))
  159.  
  160. (defun lohalf (n)            ; return low halfword
  161.   (logand n masklo))
  162.  
  163. ;; Visible functions
  164.  
  165. ;; Arithmetic functions
  166. (defun f+ (a1 a2)
  167.   "Returns the sum of two floating point numbers."
  168.   (let ((f1 (fmax a1 a2))
  169.     (f2 (fmin a1 a2)))
  170.     (if (same-sign a1 a2)
  171.     (setq f1 (fashr f1)        ; shift right to avoid overflow
  172.           f2 (fashr f2)))
  173.     (normalize
  174.      (cons (+ (car f1) (ash (car f2) (- (cdr f2) (cdr f1))))
  175.        (cdr f1)))))
  176.  
  177. (defun f- (a1 &optional a2)        ; unary or binary minus
  178.   "Returns the difference of two floating point numbers."
  179.   (if a2
  180.       (f+ a1 (f- a2))
  181.     (normalize (cons (- (car a1)) (cdr a1)))))
  182.  
  183. (defun f* (a1 a2)            ; multiply in halfword chunks
  184.   "Returns the product of two floating point numbers."
  185.   (let* ((i1 (car (fabs a1)))
  186.      (i2 (car (fabs a2)))
  187.      (sign (not (same-sign a1 a2)))
  188.      (prodlo (+ (hihalf (* (lohalf i1) (lohalf i2)))
  189.             (lohalf (* (hihalf i1) (lohalf i2)))
  190.             (lohalf (* (lohalf i1) (hihalf i2)))))
  191.      (prodhi (+ (* (hihalf i1) (hihalf i2))
  192.             (hihalf (* (hihalf i1) (lohalf i2)))
  193.             (hihalf (* (lohalf i1) (hihalf i2)))
  194.             (hihalf prodlo))))
  195.     (if (> (lohalf prodlo) round-limit)
  196.     (setq prodhi (1+ prodhi)))    ; round off truncated bits
  197.     (normalize
  198.      (cons (if sign (- prodhi) prodhi)
  199.        (+ (cdr (fabs a1)) (cdr (fabs a2)) mantissa-bits)))))
  200.  
  201. (defun f/ (a1 a2)            ; SLOW subtract-and-shift algorithm
  202.   "Returns the quotient of two floating point numbers."
  203.   (if (zerop (car a2))            ; if divide by 0
  204.       (signal 'arith-error (list "attempt to divide by zero" a1 a2))
  205.     (let ((bits (1- maxbit))
  206.       (quotient 0) 
  207.       (dividend (car (fabs a1)))
  208.       (divisor (car (fabs a2)))
  209.       (sign (not (same-sign a1 a2))))
  210.       (while (natnump bits)
  211.     (if (< (- dividend divisor) 0)
  212.         (setq quotient (ash quotient 1))
  213.       (setq quotient (1+ (ash quotient 1))
  214.         dividend (- dividend divisor)))
  215.     (setq dividend (ash dividend 1)
  216.           bits (1- bits)))
  217.       (normalize
  218.        (cons (if sign (- quotient) quotient)
  219.          (- (cdr (fabs a1)) (cdr (fabs a2)) (1- maxbit)))))))
  220.   
  221. (defun f% (a1 a2)
  222.   "Returns the remainder of first floating point number divided by second."
  223.   (f- a1 (f* (ftrunc (f/ a1 a2)) a2)))
  224.       
  225.  
  226. ;; Comparison functions
  227. (defun f= (a1 a2)
  228.   "Returns t if two floating point numbers are equal, nil otherwise."
  229.   (equal a1 a2))
  230.  
  231. (defun f> (a1 a2)
  232.   "Returns t if first floating point number is greater than second,
  233. nil otherwise."
  234.   (cond ((and (natnump (car a1)) (< (car a2) 0)) 
  235.      t)                ; a1 nonnegative, a2 negative
  236.     ((and (> (car a1) 0) (<= (car a2) 0))
  237.      t)                ; a1 positive, a2 nonpositive
  238.     ((and (<= (car a1) 0) (natnump (car a2)))
  239.      nil)                ; a1 nonpos, a2 nonneg
  240.     ((/= (cdr a1) (cdr a2))        ; same signs.  exponents differ
  241.      (> (cdr a1) (cdr a2)))        ; compare the mantissas.
  242.     (t
  243.      (> (car a1) (car a2)))))    ; same exponents.
  244.  
  245. (defun f>= (a1 a2)
  246.   "Returns t if first floating point number is greater than or equal to 
  247. second, nil otherwise."
  248.   (or (f> a1 a2) (f= a1 a2)))
  249.  
  250. (defun f< (a1 a2)
  251.   "Returns t if first floating point number is less than second,
  252. nil otherwise."
  253.   (not (f>= a1 a2)))
  254.  
  255. (defun f<= (a1 a2)
  256.   "Returns t if first floating point number is less than or equal to
  257. second, nil otherwise."
  258.   (not (f> a1 a2)))
  259.  
  260. (defun f/= (a1 a2)
  261.   "Returns t if first floating point number is not equal to second,
  262. nil otherwise."
  263.   (not (f= a1 a2)))
  264.  
  265. (defun fmin (a1 a2)
  266.   "Returns the minimum of two floating point numbers."
  267.   (if (f< a1 a2) a1 a2))
  268.  
  269. (defun fmax (a1 a2)
  270.   "Returns the maximum of two floating point numbers."
  271.   (if (f> a1 a2) a1 a2))
  272.       
  273. (defun fzerop (fnum)
  274.   "Returns t if the floating point number is zero, nil otherwise."
  275.   (= (car fnum) 0))
  276.  
  277. (defun floatp (fnum)
  278.   "Returns t if the arg is a floating point number, nil otherwise."
  279.   (and (consp fnum) (integerp (car fnum)) (integerp (cdr fnum))))
  280.  
  281. ;; Conversion routines
  282. (defun f (int)
  283.   "Convert the integer argument to floating point, like a C cast operator."
  284.   (normalize (cons int '0)))
  285.  
  286. (defun int-to-hex-string (int)
  287.   "Convert the integer argument to a C-style hexadecimal string."
  288.   (let ((shiftval -20)
  289.     (str "0x")
  290.     (hex-chars "0123456789ABCDEF"))
  291.     (while (<= shiftval 0)
  292.       (setq str (concat str (char-to-string 
  293.             (aref hex-chars
  294.                   (logand (lsh int shiftval) 15))))
  295.         shiftval (+ shiftval 4)))
  296.     str))
  297.  
  298. (defun ftrunc (fnum)            ; truncate fractional part
  299.   "Truncate the fractional part of a floating point number."
  300.   (cond ((natnump (cdr fnum))        ; it's all integer, return number as is
  301.      fnum)
  302.     ((<= (cdr fnum) (- maxbit))    ; it's all fractional, return 0
  303.      '(0 . 1))
  304.     (t                ; otherwise mask out fractional bits
  305.      (let ((mant (car fnum)) (exp (cdr fnum)))
  306.        (normalize 
  307.         (cons (if (natnump mant)    ; if negative, use absolute value
  308.               (ash (ash mant exp) (- exp))
  309.             (- (ash (ash (- mant) exp) (- exp))))
  310.           exp))))))
  311.  
  312. (defun fint (fnum)            ; truncate and convert to integer
  313.   "Convert the floating point number to integer, with truncation, 
  314. like a C cast operator."
  315.   (let* ((tf (ftrunc fnum)) (tint (car tf)) (texp (cdr tf)))
  316.     (cond ((>= texp mantissa-bits)    ; too high, return "maxint"
  317.        mantissa-maxval)
  318.       ((<= texp (- mantissa-bits))    ; too low, return "minint"
  319.        mantissa-minval)
  320.       (t                ; in range
  321.        (ash tint texp)))))        ; shift so that exponent is 0
  322.  
  323. (defun float-to-string (fnum &optional sci)
  324.   "Convert the floating point number to a decimal string.
  325. Optional second argument non-nil means use scientific notation."
  326.   (let* ((value (fabs fnum)) (sign (< (car fnum) 0))
  327.      (power 0) (result 0) (str "") 
  328.      (temp 0) (pow10 _f1))
  329.  
  330.     (if (f= fnum _f0)
  331.     "0"
  332.       (if (f>= value _f1)            ; find largest power of 10 <= value
  333.       (progn                ; value >= 1, power is positive
  334.         (while (f<= (setq temp (f* pow10 highest-power-of-10)) value)
  335.           (setq pow10 temp
  336.             power (+ power decimal-digits)))
  337.         (while (f<= (setq temp (f* pow10 _f10)) value)
  338.           (setq pow10 temp
  339.             power (1+ power))))
  340.     (progn                ; value < 1, power is negative
  341.       (while (f> (setq temp (f/ pow10 highest-power-of-10)) value)
  342.         (setq pow10 temp
  343.           power (- power decimal-digits)))
  344.       (while (f> pow10 value)
  345.         (setq pow10 (f/ pow10 _f10)
  346.           power (1- power)))))
  347.                       ; get value in range 100000 to 999999
  348.       (setq value (f* (f/ value pow10) all-decimal-digs-minval)
  349.         result (ftrunc value))
  350.       (let (int)
  351.     (if (f> (f- value result) _f1/2)    ; round up if remainder > 0.5
  352.         (setq int (1+ (fint result)))
  353.       (setq int (fint result)))
  354.     (setq str (int-to-string int))
  355.     (if (>= int 1000000)
  356.         (setq power (1+ power))))
  357.  
  358.       (if sci                ; scientific notation
  359.       (setq str (concat (substring str 0 1) "." (substring str 1)
  360.                 "E" (int-to-string power)))
  361.  
  362.                       ; regular decimal string
  363.     (cond ((>= power (1- decimal-digits))
  364.                       ; large power, append zeroes
  365.            (let ((zeroes (- power decimal-digits)))
  366.          (while (natnump zeroes)
  367.            (setq str (concat str "0")
  368.              zeroes (1- zeroes)))))
  369.  
  370.                       ; negative power, prepend decimal
  371.           ((< power 0)        ; point and zeroes
  372.            (let ((zeroes (- (- power) 2)))
  373.          (while (natnump zeroes)
  374.            (setq str (concat "0" str)
  375.              zeroes (1- zeroes)))
  376.          (setq str (concat "0." str))))
  377.  
  378.           (t                ; in range, insert decimal point
  379.            (setq str (concat
  380.               (substring str 0 (1+ power))
  381.               "."
  382.               (substring str (1+ power)))))))
  383.  
  384.       (if sign                ; if negative, prepend minus sign
  385.       (concat "-" str)
  386.     str))))
  387.  
  388.     
  389. ;; string to float conversion.
  390. ;; accepts scientific notation, but ignores anything after the first two
  391. ;; digits of the exponent.
  392. (defun string-to-float (str)
  393.   "Convert the string to a floating point number.
  394. Accepts a decimal string in scientific notation,  with exponent preceded
  395. by either E or e.  Only the six most significant digits of the integer
  396. and fractional parts are used; only the first two digits of the exponent
  397. are used.  Negative signs preceding both the decimal number and the exponent
  398. are recognized."
  399.  
  400.   (if (string-match floating-point-regexp str 0)
  401.       (let (power)
  402.     (f*
  403.      ; calculate the mantissa
  404.      (let* ((int-subst (extract-match str 2))
  405.         (fract-subst (extract-match str 4))
  406.         (digit-string (concat int-subst fract-subst))
  407.         (mant-sign (equal (extract-match str 1) "-"))
  408.         (leading-0s 0) (round-up nil))
  409.  
  410.        ; get rid of leading 0's
  411.        (setq power (- (length int-subst) decimal-digits))
  412.        (while (and (< leading-0s (length digit-string))
  413.                (= (aref digit-string leading-0s) ?0))
  414.          (setq leading-0s (1+ leading-0s)))
  415.        (setq power (- power leading-0s)
  416.          digit-string (substring digit-string leading-0s))
  417.        
  418.        ; if more than 6 digits, round off
  419.        (if (> (length digit-string) decimal-digits)
  420.            (setq round-up (>= (aref digit-string decimal-digits) ?5)
  421.              digit-string (substring digit-string 0 decimal-digits))
  422.          (setq power (+ power (- decimal-digits (length digit-string)))))
  423.  
  424.        ; round up and add minus sign, if necessary
  425.        (f (* (+ (string-to-int digit-string)
  426.             (if round-up 1 0))
  427.          (if mant-sign -1 1))))
  428.        
  429.      ; calculate the exponent (power of ten)
  430.      (let* ((expt-subst (extract-match str 9))
  431.         (expt-sign (equal (extract-match str 8) "-"))
  432.         (expt 0) (chunks 0) (tens 0) (exponent _f1)
  433.         (func 'f*))
  434.  
  435.        (setq expt (+ (* (string-to-int
  436.                  (substring expt-subst 0
  437.                     (min expt-digits (length expt-subst))))
  438.                 (if expt-sign -1 1))
  439.              power))
  440.        (if (< expt 0)        ; if power of 10 negative
  441.            (setq expt (- expt)    ; take abs val of exponent
  442.              func 'f/))        ; and set up to divide, not multiply
  443.  
  444.        (setq chunks (/ expt decimal-digits)
  445.          tens (% expt decimal-digits))
  446.        ; divide or multiply by "chunks" of 10**6
  447.        (while (> chunks 0)    
  448.          (setq exponent (funcall func exponent highest-power-of-10)
  449.            chunks (1- chunks)))
  450.        ; divide or multiply by remaining power of ten
  451.        (funcall func exponent (aref powers-of-10 tens)))))
  452.           
  453.     _f0))                ; if invalid, return 0
  454.  
  455. (provide 'float)
  456.  
  457. ;;; float.el ends here
  458.